rm(list = ls())

## packages 

library(tidyverse)
library(fixest)

## Load data + define density indicators 

dat <- readRDS('data/civey_clean.rds') %>%
  mutate(area_zip_km = area_zip / 1000000,
         pop_zip = pop_zip / 1000) %>% 
  mutate(density = pop_zip / area_zip_km,
         density_bin = ifelse(density > median(density, na.rm = T), 1, 0)) %>%
  mutate(treated_post_placebo_density = density*summer_bin,
         treated_post_placebo_density_bin = density_bin*summer_bin)

## Columns 1-4 
## treated_post_placebo_density = Pop. density x Pool season
## treated_post_placebo_density_bin Pop. density (binary) x Pool season

m1 <- feols(vote_afd ~ sw(treated_post_placebo_density, treated_post_placebo_density_bin) | user_id + date, 
            data = dat,
            split = ~ foreign_non_eu_share_bin, 
            cluster = ~ user_id) 

etable(m1)

## Columns 5-8 
## treated_post = Outdoor pool x Pool season 
## treated_post_placebo_density = Pop. density x Pool season
## treated_post_placebo_density_bin Pop. density (binary) x Pool season

m2 <- feols(vote_afd ~ treated_post + sw(treated_post_placebo_density, treated_post_placebo_density_bin) | user_id + date, 
            data = dat,
            split = ~ foreign_non_eu_share_bin, 
            cluster = ~ user_id) 

etable(m2)


